Mon. Not. R. Astron. Soc. 000, [THID (2010) Printed 24 February 2010 (MN M^X style file v2.2) 



On the angular momentum transport due to vertical 
convection in accretion discs 



Geoffroy Lesur and Gordon I. Ogilvie* 

Department of Applied Mathematics and Theoretical Physics, 
University of Cambridge, Wilberforce Road, Cambridge CBS 0WA, 



UK 



Accepted 23/02/2010. Received 



-; in original form 



ABSTRACT 

The mechanism of angular momentum transport in accretion discs has long been de- 
bated. Although the magnetorotational instability appears to be a promising process, 
poorly ionized regions of accretion discs may not undergo this instability. In this let- 
ter, we revisit the possibility of transporting angular momentum by turbulent thermal 
convection. Using high-resolution spectral methods, we show that strongly turbulent 
convection can drive outward angular momentum transport at a rate that is, under 
certain conditions, compatible with observations of discs. We find however that the 
angular momentum transport is always much weaker than the vertical heat transport. 
These results indicate that convection might be another way to explain global disc 
evolution, provided that a sufficiently unstable vertical temperature profile can be 
maintained. 
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1 INTRODUCTION 



■ As s hown by the early mode l s of | Shakura fc Sunvaevl l| 19731 ) 
and ILvnden-Bell fc Pringld (Il974l ). the structure and evo- 
lution of accretion discs are controlled by the process of 
angular momentum transport. It is often assumed that tur- 
bulence is responsible for this transport, although its na- 
ture and physical origin are barely understood. The most 
straightforward way to obtain a turbulent disc is to find a 
linear instability of the basic flow, which could drive turbu- 
lence in the nonlinear regime. However, the presence of a tur- 
bulent flow does not necessarily imply an outward transport 
of angular momentum through the disc, if the turbulence 
is sustained by sources of energy other than the Keplerian 
shear. 

The magnetorotational insta bility (MRI), poin t ed ou t 
in the context of accretion discs bv lBalbus fc Hawlevl (|l99ll ) , 
is believed to play an important role in accretion discs, 
as it produces turbulence that e fficiently tra nsports an- 
gular momentum outwards (e.g. Balbus 2003). However, 
the behaviour of the MRI in poorly ionized gas is un- 
certain. In some regions of discs, the instability might 
cease to exist because the resisti vity is too larg e, leading 
to the concept of 'dead zones' (|Gammid Il996l ). On the 
other hand, purely hydrodynamical processes might also 
be at work in discs, such as the subcritical baroclinic in- 
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stability dKlahr fc B odcnhcimcr 2003; Pet ersen et al. I l2007l ; 
ILesur fc Papaloizoul 2010t ). driven by radial entropy gradi- 
ents. 



The presence of a ver tical convective instabilit y in discs 
was initially suggested bv lLin fc Papalo izou ( 1980) as a way 
of transporting angular momentum. They pointed out that 
the temperature-dependence of the opacity in protostellar 
(or protoplanetary) discs naturally leads to convective in- 
stability in certain regimes. H owever, subsequent stud ies, 
using either linear ap proaches dRvu fc Goodman! Il992h o r 
nonlinear simulations (|Cabotl [l996l ; Stone fc Balbud 19961 ). 
predicted a weak inward transport associated with convec- 
tion, leading to the conclusion that this process was irrele- 
vant for the problem of accretion disc dynamics. Note that 
these simulations were performed in a 'weakly nonlinear' 
regime with relatively low Reynolds and Rayleigh numbers 
l|Cabotll 19961 ). 



In this letter, we challenge these conclusions using high- 
resolution spectral methods and a simplified setup for the 
convection problem. In we introduce our model, the rel- 
evant dimensionless parameters and the numerical setup. 
We present simulations of convection with free-slip vertical 
boundary conditions in and of homogeneous convection 
in f|4] Finally, we briefly discuss these results and their im- 
plications in jj5] 
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2 MODEL AND EQUATIONS 
2.1 Model 

In the following, we consider a model flow in which ro- 
tation, shear and thermally induced convection are in- 
cluded. Instead of considering the full disc problem, with 
its particular density and tem perature structure, we use 
the shearing box approximation (lHawlev et a l. 1995; Balbusi 
120031 ; iReeev fc Umurhar]|200Sl ). As a simplification, we as- 
sume the flow is incompressible and we introduce ver- 
tical stratification with in the Boussinesq approximation 
l|Spiegel fc Veroiiisl [i960). Although this approximation is 
not fully justified in a real disc, it allows us to explore the 
regime of strongly turbulent sheared convection in a more 
controlled way. In this letter we explore the way turbulent 
convection transports angular momentum, but we do not 
include the generation of an unstable temperature profile 
due to viscous heating. Instead, we impose a 'background' 
stratification profile as prescribed by the Boussinesq approx- 
imation. Therefore, one has to keep in mind that convection 
is not generated self-consistently in our model. 

The shearing-box equations are found by considering a 
Cartesian box centred at r = Ro, rotating with the disc at 
angula r velocity O = Q,(Rp) . Defining r—Ro — > x and Ro(j> — } 
y as in lHawlev et alj l|l995h . one obtains the following set of 
equations: 

d t u + u ■ Vu = - VII - 2f2 X u 

+2Q.Sxe al - AN 2 8e z + vAu, (1) 
d t 8 + u-V8 = u z /A + X A8, (2) 
V-u = 0, (3) 

where u is the total velocity, 8 = 5p/p is the perturbation 
of the density logarithm (or entropy, as pressure perturba- 
tions are much smaller in the Boussinesq approximation), 
v is the kinematic viscosity and x i s the thermal diffusiv- 
ity. In these equations, we have defined the mean shear 
S = —rd r Q, which is set to S = (3/2)Q for a Keplerian disc. 
The generalised pressure II is calculated by solving a Pois- 
son equation derived from the incompressibility condition. 
For dimensional consistency with the traditional Boussi- 
nesq approach, we have introduced a stratification length 
A = —g/N 2 , where g is the vertical component of the grav- 
ity. Note, however, that A disappears from the dynamical 
properties of this system of equations as one can renormal- 
ize the variables by defining 8' = A8. The stratification itself 
is controlled by the Brunt -Vaisala frequency TV, defined for 
a perfect gas by 

iV 2 = -^f 1-1* (A (4) 

7/5 dz dz \pi j w 

where P and p are the pressure and density of the back- 
ground equilibrium profile and 7 is the adiabatic index. 
Physically, TV corresponds to the oscillation frequency of a 
fluid particle displaced vertically from its equilibrium po- 
sition. When the flow is convectively unstable, one has 
TV 2 < 0. In our model, we assume TV is constant, which 
formally corresponds to a local model in z. 

One can easily check that the velocity field U = —Sxe y 
is a steady solution of equations H])-©. In the following we 
consider the evolution of the perturbations v (not necessarily 
small) of this profile defined by v = u — U. 



In the following, we us e shearing-sheet boundary con- 
ditions l|Hawlev et al.|[l995l ) in the x direction and periodic 
boundary conditions in y. In the z direction, we use either 
free-slip boundary conditions (Jj3j, imposing v z = 8 — and 
d z v x = d z Vy — 0, or periodic boundary conditions (Q. The 
free-slip vertical boundary conditions correspond to a rigid 
wall with no tangential stress and a fixed temperature. Sim- 
ilar boundary conditions for the velocity field were used by 
ICabotJ (| 19961 ). 

2.2 Dimensionless numbers 

To simplify the analysis of convection in sheared flow, we 
will use three dimensionless numbers defined as follows: 

• The Richardson number Ri = —N 2 /S 2 compares buoy- 
ancy forces to the mean shear. In a real disc, this number 
is determined by the balance between viscous heating and 
vertical heat transport. Assuming a typical velocity v due to 
convective motions, we define the turbulent viscosity using 
mixing length theory vt ~ SvH , where <5 is an efficiency fac- 
tor and H is the disc thickness. The turbulent heat transport 
is on the other hand F ~ pv 3 . Assuming thermal balance, 
we get pv 3, j H ~ pu t 0. 2 leading to v 2 ~ Sc 2 where c s is the 
local sound speed. Since turbulence is driven by the unstable 
entropy gradient, we have in first approximation v ~ \N\H, 
St.) that Ri ~ S 1/2 . As we will see in the following, we gener- 
ally get S <C 1, meaning that the regime Ri < 1 is relevant 
to accretion discs. 

• The Prandtl number Pr = v/x compares viscous dis- 
sipation to thermal diffusion. Using standard values for the 
viscosity and thermal conductivity of H2, relevant to a pro- 
toplanetary disc, one has 

where re is the opacity. 

• The Rayleigh number Ra = — N 2 L 4 /xi / compares 
buoyancy forces to dissipation processes. In a disc, one gets 




In the following, we will assume Pr — 1 for simplicity. A 
more complete study including the effect of varying Pr will 
be the subject of a subsequent paper. The length L intro- 
duced in these definitions is assumed to be a typical di- 
mension of the box and should be of the order of the disc 
thickness. In the following, we define L as the smallest box 
length with L = L z in i|3]and L — L x in |J4] When not explic- 
itly mentioned, the unit of time is the shear timescale S _1 
which is equal to (37r) -1 orbital periods. In the following, 
all of our simulations are computed in the shear dominated 
regime (Ri < 1) which is expected in discs. 

2.3 Turbulent transport and energy budget 

By analogy with IStone fc Balbusi l|l996l ) , we define the di- 
mensionless transport of angular momentum through the a- 
like coefficient a = (v x v y ) / S 2 L 2 , (■) being a time and space 



Vertical convection in discs 3 




100 200 300 400 500 

t 



Figure 1. Turbulent transport as a function of time for run Wl 
(lower curve) and W8 (upper curve). The modification of the 
average turbulent transport is due to the Rayleigh number (see 
text). 

average. Since S = 1 and L = 1 in all the simulations pre- 
sented below, we will simply write a — (v x v y ) in our units. 

Multiplying the y component of JTJ by v y , integrating 
over the volume of the box and using the boundary condi- 
tions presented above, we obtain the energy budget of the 
fluctuations in the y direction, 

d t vl — UdyVy + (S - 2Q)WU^ - v\Vv y \ 2 , (5) 

where X denotes the volume average of X. Considering a 
quasi-steady turbulence and assum ing the pressure-strain 
correlation Tld y v y to be negligible, IStone fc Balbusl (Il996l ) 
concluded from this equation that (S — 2Q.)v x v y > 0, lead- 
ing to a < in a disc. This hypothesis was motivated by 
numerical simulations (|Stone fc Balbusl 1 19961 ; ICabotl 1 19961 ) 
in which the turbulence was weakly non-axisymmetric, lead- 
ing to the conclusion that y derivatives should be negligible. 
It is however well known that the pressure-strain correla- 
tion tensor is responsible for the redist ribution of energy in 
fully developed anisotropic turbulence (Spczialc 1991; Pope 
2000). Therefore, when the flow becomes strongly nonlinear 
(i.e. at large Ra), there is no reason to assume that this term 
is negligible. 

2.4 Numerical model 

We use the Snoopy code to solve equations (HJ) — dSJl . This 
code uses a 3D Fourier representation of the flow with a 
third-order timestepping scheme to have an excellent accu- 
racy at all scales. It has been used in sev eral studies of hydro- 
dynamic an d MHD turbulence (see e.g. iLesur fc Longarettil 
120051 . 120071 1. When working with free-slip boundary condi- 
tions, a sine-cosine decomposition is used instead of the 
Four ier basis in the z direction (see e.g. Ca ttaneo et al.l 
l2003h . 



3 CONVECTION WITH WALLS 

In the first set of simulations, we introduce free-slip bound- 
ary conditions in the z direction. These boundary conditions 
prevent the development of large-scale vertical flows which 
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Figure 2. Correlation between the turbulent transport and Ra 
for Ri = 0.4 (plain line) and Ri = 0.2 (dashed line). 

are unrealistic for a disc (see In the particular setup con- 
sidered here, convection starts at Ra ~ 4100 for Ri = 0.4 
and Ra ~ 6900 for Ri = 0.2. 

In order to study the dependence of the turbulent trans- 
port of angular momentum on our dimensionless parame- 
ters, we have carried out several simulations which are pre- 
sented in Tab.[T] Average values are computed from simula- 
tions run for 500 <S _1 , excluding the first 100 5 _1 from the 
averaging procedure to eliminate the influence of the initial 
conditions. As shown in Fig. [1] the simulations are suffi- 
ciently homogeneous on these timescales for this average to 
be meaningful. 

3.1 Rayleigh number and transport 

The most striking observation in these simulations is the cor- 
relation between Ra and the turbulent tra nsport a. In par- 
ticula r, for Ra < 10 6 we recover the result of lStone fc Balbusl 
l|l996h of a small and negative transport (runs W1-W4 and 
W9-W12). However, for sufficiently large Ra, the sign of the 
transport is reversed. To demonstrate this effect, we show 
in Fig. [2] the dependence of the turbulent transport on Ra, 
for Ri = 0.4 and Ri = 0.2. 

It is possible to deduce a general scaling for the trans- 
port from these results: 

a ~ -1.8 x 10~ 3 + 2.9 x 10~ 4 \og 10 (Ra) (6) 
for Ri = 0.4 and 

a ~ -5.4 x 10" 4 + 7.2 x 10~ 5 log 10 (i?a) (7) 

for Ri — 0.2. Surprisingly, these scalings predict a transport 
of the order of 10 -3 for values of Ra which could be expected 
in a disc. However, this conclusion is challenged by the fact 
that at large Ra, the dependence of a on Ra should become 
weaker since the length-scales associated with transport are 
then much larger than the diffusive scales and the diffusion 
coefficients may be expected to become less important. A 
weak effect of this kind is already observed in Fig. [5] for 
Ra > 10 7 . 

Interestingly, changing the Rayleigh number does not 
significantly modify the turbulent kinetic energy of the sys- 
tem. This clearly demonstrates that the correlation observed 
above is not due to a modification of the turbulent activity in 
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Run 


Resolution 


Box size 


no 


hi 




„. \ / c2 r 2 

(VxVy)/b = 


a 


A{6v z )/SL 2 




Tlx X Thy X Tlz 




















Wl 


256 x 256 x 128 


2x4x1 


1 X 10 5 


0.4 


5.7 x 10~ 


3 


-3.9 x 10- 


-4 


7.7 x 10- 


3 


W2 


256 x 256 x 128 


2x4x1 


3 x 10 5 


0.4 


7.0 x 10" 


3 


-2.3 x 10- 


-4 


7.4 x 10- 


3 


W3 


256 x 256 x 128 


2x4x1 


6 x 10 5 


0.4 


6.9 x 10" 


3 


-1.4 x 10- 


-4 


7.0 x 10" 


3 


W4 


256 x 256 x 128 


2x4x1 


1.2 x 10 6 


0.4 


7.2 x 10" 


3 


-4.6 x 10" 


-5 


6.2 x 10" 


3 


W5 


256 x 256 x 128 


2x4x1 


2.5 x 10 6 


0.4 


7.2 x 10" 


3 


5.8 x 10" 


-5 


5.7 x 10" 


3 


W6 


256 x 256 x 128 


2x4x1 


5 x 10 6 


0.4 


7.1 x 10" 


3 


1.4 x 10- 


-4 


5.1 x 10" 


3 


W7 


256 x 256 x 128 


2x4x1 


1 x 10 7 


0.4 


6.7 x 10- 


3 


2.1 x 10- 


-4 


4.6 x 10- 


3 


W8 


256 x 256 x 128 


2x4x1 


1.4 x 10 7 


0.4 


6.4 x 10" 


3 


2.4 x 10- 


-4 


4.3 x 10- 


3 


W9 


256 x 256 x 128 


2x4x1 


3 x 10 5 


0.2 


2.1 x 10- 


3 


-1.6 x 10- 


4 


3.9 x 10- 


3 


W10 


256 x 256 x 128 


2x4x1 


1.2 x 10 6 


0.2 


2.2 x 10" 


3 


-1.0 x 10- 


-4 


3.7 x 10- 


3 


Wll 


256 x 256 x 128 


2x4x1 


5 x 10 6 


0.2 


2.2 x 10- 


3 


-5.3 x 10- 


-5 


3.1 x 10- 


3 


W12 


256 x 256 x 128 


2x4x1 


2 x 10 7 


0.2 


2.0 x 10- 


3 


-8.7 x 10- 


-6 


2.5 x 10- 


3 


W13 


512 x 512 x 256 


2x4x1 


8 x 10 7 


0.2 


1.8 x 10- 


3 


3.0 x 10- 


-5 


1.9 x 10- 


3 


W14 


512 x 512 x 256 


2x4x1 


3.2 x 10 s 


0.2 


1.7 x 10" 


3 


5.8 x 10- 


-5 


1.6 x 10- 


3 



Table 1. List of simulations with free-slip boundary conditions. The turbulent kinetic energy (6 th column), angular momentum transport 
(7 th column) and heat transport (8 th column) are computed from simulations. 
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Figure 3. Energy budget in the y direction (eq. [5} from run W7. 

the flow but results instead from a change in the anisotropy 
of the turbulence. Note also that the turbulent heat trans- 
port (8v z ) decreases as Ra increases. However, we always 
find (0Vz) S> (v x v y ), showing that convection is inefficient 
at transporting angular momentum. This point will be ex- 
amined further in iJS] 

3.2 Energy budget 

To check if the pressure-strain correlation is important in 
our simulations, we plot in Fig. [3] the time history of each 
term of equation (JSJ for run W7. As shown in this plot, the 
energy balance at high Ra is dominated by an equilibrium 
between the pressure-strain correlation and the dissipation 
term. Moreover, the turbulent transport term is negligible in 
this situation. Thi s clearly indicate s that at high enough Ra, 
the hypothesis of IStone fc Balbusl l|l99fj ) is no longer valid 
and the turbulent transport of angular momentum can be 
directed outward. 

The presence of a large pressure-strain correlation term 
in equation ((SJ) shows that non-axisymmetric structures are 
playing an important role at large Ra. To isolate the role 
played by these non-axisymmetric structures, we plot in 
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Figure 4. ID spectra of the pressure-strain correlation from run 
Wl (plain line), W5 (dashed line) and W7 (dash-dotted line). 

fig.|4]the ID spectrum of the pressure-strain correlation. As 
expected, we find no contribution from axisymmetric modes 
(k y = 0). However, as Ra becomes larger, the contribution 
from smaller-scale non-axisymmetric modes becomes more 
important. This effect is evident in run W7 where contribu- 
tions from scales as small as L y /20 are still important. 

This result clearly indicates that as the flow becomes 
turbulent, it develops small-scale motions which tend to re- 
distribute the energy fluctuations in the three directions. 
This is consistent with the picture of anisotropic turbulence 
becoming isotropic at small scales thanks to the pressure- 
strain correlation. It also indicates that the results o f lCabotl 
jl996t ) and lStone fc Balbusl (|l996t ) were due to a too small 
Ra, an issue already mentioned bv lCabotl (|l996t ). 



4 HOMOGENEOUS CONVECTION 

To test the possible limitations of the free-slip boundary 
conditions used above, we have also considered the prob- 
lem of homogeneous convection, using periodic boundary 
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Run 


Resolution 


Box size 


Ra 


Ri 


{\v 2 )/S 2 L 2 


(VxVy)/S 2 L 2 = 


a 


A{8v z )/SL 2 




n x x n y X n z 


Lx X Liy X Lz 














HI 


64 x 64 x 128 


1x1x2 


10 4 


0.1 


2.8 x 10" 2 


-3.4 x 10" 


-4 


1.1 x lO" 1 


H2 


64 x 64 x 128 


1x1x2 


4 x 10 4 


0.1 


1.8 x 10" 2 


9.6 x 10" 


-4 


6.2 x 10~ 2 


H3 


128 x 128 x 256 


1x1x2 


4 x 10 5 


0.1 


1.6 x 10" 2 


9.0 x 10- 


-4 


4.2 x 10~ 2 


HI 


128 x 128 x 256 


1x1x2 


4 x 10 6 


0.1 


1.9 x 10" 2 


1.6 x 10" 


-3 


3.7 x 10~ 2 


H5 


128 x 256 x 256 


1x2x2 


4 x 10 5 


0.1 


1.5 x 10" 2 


5.3 x 10" 


-4 


4.0 x 10~ 2 


H6 


128 x 256 x 256 


1x2x2 


4 x 10 s 


0.1 


1.3 x 10" 2 


9.0 x 10 


-4 


2.8 x 10~ 2 



Table 2. List of simulations with periodic boundary conditions in the z direction. 



conditions in the vert ical direction applied to v and 6 
jBorue fc Orszad 1 1997ft . In this context, the flow becomes 
convectively unstable for Ra > 1558, developing first "eleva- 
tor" modes with a purely vertical motion u 2 oc cos(kx% + <f>)- 
Since these modes are also exact nonlinear solutions of equa- 
tions CO)-©, they can grow exponentially for long times, 
leading to "spikes" in the turbulent energy and transport. 
These exponential growth phases are unrealistic, as they 
should be broken by secondary instabilities or the effects of 
density stratification, which are neglected in the Boussinesq 
approximation. A way to limit the formation and amplifica- 
tion of these modes is to consider "tall" boxes with L z > L x , 
so that secondary instabilities are more efficient at break- 
ing up the elevator modes. Note that the phenomenology of 
these elevator modes is very similar to that of channel mode 
solutions fou nd in the context of th e magnetorotational in- 
stability (see lGoodman fc Xu|[l99i '). 

We have used this tall box setup with several values 
of Ra and L y in the regime of weak convection, Ri = 0.1 
(Tab. [2]). Looking at the time series (Fig. [5]), we note the 
presence of large peaks of transport at small Ra, with both 
positive and negative values. These peaks are associated 
with the growth of elevator modes to high amplitude and 
their break-up through secondary instabilities. We note how- 
ever that this behaviour disappears at larger Ra, although 
large fluctuations remain, leading to a positive-definite tur- 
bulent transport with an average value a ~ 10~ 3 . 

Modifying the size of the box in the y direction does 
not change dramatically the results, and on average, we find 
a transport that is stronger in the homogeneous case than 
with free-slip boundary conditions. Note also that the aver- 
age turbulent kinetic energy is larger in homogeneous runs 
(Tab. [3J) which partly explains the higher turbulent trans- 
port observed in those cases. 



5 DISCUSSION 

In this letter, we have shown that turbulent convection, if 
present in accretion discs, could lead to an outward trans- 
port of angular momentum, provided that the Rayleigh 
number is large enough (i.e. molecular viscosity and ther- 
mal diffusivity are small enou gh). We also conjectured that 
the p revious negative results jCabod 119961 ; IStone fc Balbusl 
Il996ft were due to a too small Rayleigh number, leading to 
nearly axisymmetric turbulence which cannot transport an- 
gular momentum outward. 

Despite this positive result, this letter does not demon- 
strate that convection is actually present in discs and is re- 
sponsible for the observed accretion rates. To demonstrate 




Figure 5. Turbulent transport as a function of time for run H2 
(dashed line) and H4 (plain line). 

this point, one probably has to include viscous heating and 
a realistic model of disc stratification, which are beyond the 
scope of this letter. 

We note that the use of the Boussinesq approximation 
with constant stratification is a strong simplification of the 
problem. However, we always find v < 0.1SL, showing that 
our incompressible approximation is valid. Moreover, we be- 
lieve that the Rayleigh number effect found in this work will 
remain valid with a more realistic vertical profile, as small 
scales are weakly influenced by the global (vertical) struc- 
ture. 

Interestingly, we find that the turbulent transport of 
angular momentum is always much weaker than the turbu- 
lent heat transport (typically (9v z ) ~ 20a). This is to be 
expected as turbulent convection is not shear driven, but 
it also indicates that convection is inefficient at transport- 
ing angular momentum. This leads us to conjecture that it 
might not be possible to generate vertical convection self- 
consistently with viscous heating only. Instead, one might 
have to consider an additional heating source which could 
be due to radiation or chemical reactions. 
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